Example of 20 hypothetical 1ha restoration plots (10,000m2, 100m*100m) distributed at random at three reefs in the Townsville region (Davies Reef, Big Broadhurst Reef, Little Broadhurst Reef). See code and “Visualise restoration plots” for workflow.

# data packages
library(tidyverse)
library(units)
library(foreach)

# spatial packages
library(sf)

# mapping packages
library(tmap)
library(leaflet)

# create function to get make polygons for plots
set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {

  # Calculate the coordinates of the rectangular polygon
  x <- sf::st_coordinates(input)[1, 1]
  y <- sf::st_coordinates(input)[1, 2]

  # set parameters
  x_min <- x - (width / 2)
  x_max <- x + (width / 2)
  y_min <- y - (length / 2)
  y_max <- y + (length / 2)

  polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
    sf::st_sfc(crs = 20353)

  return(polygon)
}



habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")

cairns_geomorphic <- st_read("/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson") %>%
   sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
   dplyr::group_by(class) %>%
   dplyr::mutate(habitat_id = paste0(
   gsub(" ", "_", class), "_",
   sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
   sf::st_make_valid() %>% 
#  mutate(habitat_id=as.factor(habitat_id)) %>%
  dplyr::group_by(habitat_id, class) %>% 
  summarise(geometry = st_union(geometry)) %>% 
  mutate(area = round(st_area(geometry))) %>%
  filter(class %in% habitats) %>%
  drop_na(class) 
## Reading layer `Moore Arlington' from data source `/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson' using driver `GeoJSON'
## Simple feature collection with 3768 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 145.9182 ymin: -16.89574 xmax: 146.2582 ymax: -16.61326
## Geodetic CRS:  WGS 84
# find restorable habitats
cairns_sheltered_sites <- filter(cairns_geomorphic, class %in% c("Sheltered Reef Slope", "Outer Reef Flat")) 
# identify neighbours
neighbors <- st_touches(cairns_geomorphic)
# match pairs
border_rows <- purrr::map2(neighbors, cairns_geomorphic$class, ~ {
  any(cairns_geomorphic$class[.x] == "Sheltered Reef Slope") & (.y == "Outer Reef Flat")
})

# Filter rows
indices <- which(unlist(border_rows))
result <- cairns_geomorphic[indices, ]

# combine outer and sheltered
cairns_restorable <- rbind(
  filter(cairns_geomorphic, class %in% c("Sheltered Reef Slope")),
  result
)


cairns_plots_centroids <- cairns_restorable %>%
  filter(area > set_units(10000,"m^2")) %>%
  filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>% 
  drop_na(class) %>%
  st_centroid() %>% 
  st_cast("POINT")


cairns_plot_outputs <- foreach(i=1:nrow(cairns_plots_centroids), .combine="c") %do% {
  tmp <- set_restoration_plot_centres(cairns_plots_centroids$geometry[i], width=100, length=100)
  tmp
}

Map restoration plots

Select layers using the layercontrol on the left, use [ ] for full-screen viewing.

set.seed(1)

cairns_plot_outputs <- cairns_plot_outputs |> st_as_sf() |> mutate(id=paste0("Plot_ID_",seq(1,n(),1))) 

cairns_plot_outputs2 <- cairns_plot_outputs[sample(nrow(cairns_plot_outputs), 20), ] # subset to 20


# visualise with thematic maps
map <- tm_view() +
#tm_tiles("Esri.WorldImagery", group = "[Satellite map]", alpha = 0.5) +

#---------- all habitats--------@
tm_shape(cairns_geomorphic |> mutate(class=as.factor(class)),  name = "Reef habitats", id="area") +
  tm_borders(col = "black", lwd = 0.2) +
  tm_fill("class", name = "Benthic habitats", Title = "Benthic habitats", id="area", palette = c("Plateau" = "lightskyblue1", "Back Reef Slope" = "darkcyan",
                                                                    "Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
                                                                    "Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
                                                                    "Reef Crest" = "coral3"), alpha = 0.6) +
#---------- restorable habitat--------@
tm_shape(cairns_restorable |> mutate(class=as.factor(class)), id="area", name = "Restorable habitats") +
  tm_borders(col = "white", lwd = 0) +
  tm_fill("class", title="Restorable habitats", name = "Benthic habitats", id="area", palette = c("Sheltered Reef Slope" = "darkslategrey",  "Outer Reef Flat" = "darkgoldenrod3"), alpha = 0.6) +

#---------- restoration plots--------@
tm_shape(cairns_plot_outputs2, id="plot_id", name = "Restoration plots (1ha)") +
    tm_fill(col="red", alpha=0.6) +
    tm_borders(col="black") +
#---------- scale--------@
tm_scale_bar(width=200)

  map |> tmap::tmap_leaflet() |>
    leaflet::addProviderTiles('Esri.WorldImagery', group = "Satellite map", options=leaflet::providerTileOptions(opacity=0.6, maxNativeZoom=18,maxZoom=100)) |>
#    leaflet::addProviderTiles('Esri.WorldTopoMap',  group = "<b> [Restoration]</b> habitats", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=10)) |>
    leaflet::addLayersControl(position="topleft", overlayGroups=c(
      "Satellite map", 
      "Reef habitats",
      "Restorable habitats",
      "Restoration plots (1ha)"),
                              options=leaflet::layersControlOptions(collapsed = FALSE)) |>
     leaflet::hideGroup(c("Reef habitats")) |> 
    leaflet.extras::addFullscreenControl(position = "topleft", pseudoFullscreen = FALSE)

Summary statistics

Habitat Area TwentyPlots
Total coral area 85598086 0.2%
Total restorable area 19098405 1%